The Heterogeneous Multi-Scale Method 

Weinan E* Bjorn Engquist^ 

February 2, 2008 



Abstract 

The heterogeneous multi-scale method (HMM) is a general strategy for dealing 
with problems involving multi-scales, with multi-physics, using multi-grids. It not 
only unifies several existing multi-scale methods, but also provide a methodology for 
designing new algorithms for new applications. In this paper, we review the history 
of multi-scale modeling and simulation that led to the development of HMM, the 
methodology itself together with some applications, and the mathematical theory of 
stability and accuracy. 

1 Introduction 

In the past several years, there has been an explosive growth of interest on numerical com- 
putations for problems involving multi-scales, often with multiple levels of physical models 
and use multi-grids. Applications of these ideas are found in many different areas, including 
coupling quantum mechanics with molecular dynamics |L2| [73|, [75|, coupling atomistics with 
continuum theory |)6|, [l], |55|, fx| 11, 24, 45, 37], coupling kinetic theory with continuum 



theory || ^0], [4^, [72], coupling kinetic Monte Carlo methods with continuum theory [|58 
homogenization theory |5|, |28], [5], |60], [|T|, and coarse-grained bifurcation analysis |67|, |57| . 
^From the point of view of numerical analysis, it is natural to ask whether these seemingly 
different applications can be put into a common framework, and whether a general theory 
of stability and accuracy can be provided. Such a theory should help us to improve existing 
methods, design new ones, and extend their applicability to other problems. 

Finite element method provided an example of a success story of this kind. The practice 
of finite element methods was started by structural engineers on very specific applications. 
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However, the work of mathematicians at the end of the 60's and early 70's put finite element 
method on a much more solid and general framework. This framework provided a thorough 
understanding of existing methods, suggested improvements and as well as extension to other 
areas such as fluid mechanics, electromagnetism, etc. 



In this article, we will review the work in [^T|] that aimed at providing a general framework 



as well as a stability and accuracy theory for numerical methods involving multiscales, multi- 
grids and multi-physics. For convenience and to emphasize the multi-physics nature of these 
methods, we will call them Heterogeneous multiscale methods. In contrast, typical multi- 
grid methods are homogeneous in the sense that they use the same model at different levels 
of grids and are aimed at efficiently resolving the details at the finest grids. 

To begin with, let us make some remarks about problems with multiple scales. Such 
problems are found everywhere around us. A classical example that has been extensively 
studied in the mathematics literature is the problem of homogenization. 

-V'(a^)Vtt £ (x))=/(x), xen (1) 

with Dirichlet boundary condition u E \gn = 0, for example ||. Here the multiscale nature 
is reflected in the coefficients a (x, |) , £ -C 1. In the simplest models a(x, y) is assumed to 
be periodic in y, say with period I = [0, l] d , d is the dimension of the physical space. ([[]) 
can be used for modelling transport properties in a medium with microstructure, and the 
oscillatory nature of a is used to model the microstructures in the medium. 

Traditionally problems of this type are dealt with either analytically or empirically by 
finding effective models that eliminate the small scales. For the homogenization problem 
(ED), this means replacing (ffl) by a homogenized equation [0, or effective equation 



-V • [A{x)VU{x)) = f{x) xett. (2) 

The solution U of this equation approximates the behavior of u £ averaged over length scales 
that are much larger than e but smaller than the slow variations of a and /. In the special 
case when d — 1, A(x) is given simply by 



An althernative to such analytical methods is the empirical modelling. As an example, 
let us consider simple incompressible fluids. Let u be the velocity field. Mass and momentum 
conservation gives, 

u t + (u- V)w + -Vp = V • r d , V • u = (4) 
P 

where p is the (constant) density, is the viscous stress, which is a macroscopic idealization 
of the internal friction forces due to the short-ranged molecular interactions, must be 
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modeled in order to close the equation (f|). The simplest empirical model is to assume that 
Td is linearly related to D — V "~H V ") ; the ra ^ e G f strain tensor. Using homogeneity, isotropy 
and incompressibility gives the constitutive relation 

r d = vD (5) 

where v is the viscosity of the fluid. In this model, all molecular details are lumped into 
a single number v. It is quite amazing that such a simple ansatz describes very well the 
behavior of simple liquids in almost all regimes. 

Despite all these successes, such traditional approaches also have their limitations. Al- 
though a nice set of equations can be derived rigorously for the effective coefficients A(x) 
in (H), they are not explicit and evaluating them involves a considerable amount of work. 
Finding empirical constitutive relations for complex fluids such as polymeric fluids or liquid 
crystals has also proven to be a difficult task. In general the constitutive relations tend to 
contain many empirical parameters and their accuracy is also often in serious doubt. 

Such problems are not limited to hydrodynamics where constitutive relations are needed, 
but is common to most modeling process. In molecular dynamics, one needs to model, often 
empirically, the atomic potentials. In kinetic Monte Carlo methods, one needs to model 
the transition rates. In kinetic equations, one needs to model the collision cross-section. 
In nonlinear elasticity, one needs to model the stored-energy functional. In typical mean 
field theories, one needs to model the effective local fields and the free energy functional. In 
general, the empircal models work well for relatively simple systems, but lose their accuracy 
for complex systems. 

In the last decade, a new approach, the "first principle-based" approach, has been ve- 
hemently pursued in various areas of applications. The basic idea is to replace empirical 
models by coupling with direct numerical simulations at a more microscopic level. Some of 
the best known examples of this approach include: 

1. Ab initio molecular dynamics [12|. Here one replaces the empirical atomic potential in 
molecular dynamics by explicit calculations of the electronic structures. The Car-Parrinello 
method is a practical way of implementing this idea [12 |. 

2. Quasi-continuum method ^2[. This is a method for doing nonlinear elasticity 



calculations without the need of a stored-energy functional. Instead, the stored energy is 
computed directly from atomic potentials using the Cauchy-Born rule. We will return to 
this later. 

3. The Gas-kinetic scheme [[71 . This is a method for doing gas dynamics calculations 



using directly the kinetic equations instead of the hydrodynamic equations. Since it played 
an important role in the framework developed in we briefly review the important steps 
here. 

Given {p n , u n , T n }, the density, velocity and temperature at time step n at each cell, the 
corresponding values at the next time step, {p n+1 ,u n+1 ,T n+1 } are computed by: 
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Step 1. Reconstruction. From {p n ,u n ,T n }, reconstruct f n , the one particle phase space 
distribution function near the cell boundaries. 

Step 2. Solve the kinetic equation with initial data f n near the cell boundaries. In [ 7I |, 
the kinetic equaiton is chosen as the BGK equation 

ft + (u-V)f = - e (X (p , u , T) -f) 

where X( p , u ,t) is the local Maxwellian associated with (p,u,T). 

Step 3. Compute the average density, momentum and energy fluxes at the cell boundaries, 
from which one computes {p n+1 , u n+1 , T n+1 }. 

This procedure is an illustration of several of the key ingredients that we use in the general 
framework that we will discuss below: the selection of an overall macroscale scheme which 
in the present example is the finite volume method; the estimation of the macroscale data, 
here the flux, using the Godunov procedure which consists of the steps of reconstruction, 
microscale evolution and averaging; the cost reduction at the microscale evolution step by 
restricting to a small subset of the computational domain. 

The examples we discussed so far belong to the class of problems referred to as type B 
problems in |2lJ where microscopic models are used to supply a closure to the macroscopic 
models. Another wide class of problems, called type A problems in |H] consist of problems 
with defects, interfaces or singularities, for which conventional macroscopic models are ac- 
curate enough away from the defects, and more detailed microscopic models are necessary 
near the defects. Type A problems are found in crack propagation, contact line dynamics, 
triple junctions, grain boundary motion, dislocation dynamics, etc. 

This short review of existing work is by no means comprehensive. For the convenience 
of the reader, we include at the end an extensive list of references. 



2 Relations between Macroscopic and Microscopic Mod- 
els 

Let us first fix the notations. We will denote the macroscopic and microscopic state variables 
as U and u, defined on D and V, respectively. Typically D is the physical space. As we 
will explain below using examples, it is convenient to view T> loosely as a fiber bundle over 
D, where the fiber T> x over x G D, is roughly the space of microstructures at x. The 
macroscopic and microscopic state variables are connected by a compression operator Q, 
and a reconstruction operator R 

Qu = U, RU = u, QR = I 

where I is the identity operator. Typically there is a natural way to define Q. But R is 
certainly not unique. 
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To illustrate these notions and notations, let us consider a few examples. 

1. For the first example we will consider the case when the macroscopic variables are 
the hydrodynamic variables of density, velocity and temperature (p,u,T), and the 
microscopic variable is the one particle phase space distribution function /. In this 
case D is the physical space-time domain of interest, T> = D x R 3 , T> x = R 3 , the 
momentum space which represents the microstructures in this problem. Q is defined 
by 

p(x,t) = / f(v,x,t)dv, u(x,t) = / f(v,x,t)vdv, T(x,t) — — / f(v, x, t)\v — u\ 2 dv 
Jb? Jr 3 3p J R 3 

2. For our second example, we will take the microscopic model to be the kinetic mod- 
els of rod-like molecules in liquid crystals, [|18| where the microscopic variable is the 
orientation-position distribution function f(x,m,t), (x,t) G D, m G S 2 , the unit 
sphere. Here T> = D x S 2 , T> x = S 2 . The macroscopic model is the Landau-de Gennes 
type models with tensorial order parameter S which is our macroscopic variable. Q is 
defined as 

S(x,t) = / [m®m 1 ) f(x,m,t)dm 
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3. For the third example, we take the standard homogenization problem 



U+ + ( a ( x, - ) u e ) =0 



where a(x, y) is smooth and periodic in y with period 1. In this case, the macroscopic 
variable will be the local space-time averages of u e : 



1 f 

U(x,t) = — — / u £ {x + y,t + s)dyds 
\^\ Jc 



where |C| denotes the area of c on which the averaging is taken. D x = (x, t) + C . The 
size of C should be larger than e. 

4. Our last example is front propagation described by Ginzburg-Landau equations 

u\ = Au £ + ^u £ {l-{u £ ) 2 ) 

To define the macroscopic variables, observe that u £ is close to ±1 in most of the 
physical domain D, except in a thin region of thickness 0(e) where sharp transition 
between ±1 takes place. Since the fast reaction term vanishes at three values —1, 0, +1, 
it is natural to define U by: 

1 ifu £ {x,t)>0 
U(x,t) = Qu £ (x,t) = { iiu £ (x,t) = 

-1 ifw £ (x,t)<0 
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An equivalent definition of U is via the level set of u e . 
More examples are disscussed in JZTJ. 

In the following we will concentrate on problems of type B, namely problems for which 
there exist a set of macroscopic variables that obey closed macroscopic models, but the 
macroscopic models are not explicitly available. We will describe general computational 
methodologies that enable us to do numerical computations efficiently based on the under- 
lying microscopic models. We will comment on problems of type A at the end. 



3 Abstract Formulations 

A key component of HMM is the estimation of macroscale data using the microscale model. 
To see how this can be done, it is helpful to first give an abstract formulation of the macro- 
scopic model in terms of the microscopic model. We start with variational problems. 
Consider a microscopic variational problem 

mme(u) (6) 

where u) is some function space over £>, the physical space. Let Q be the compression 
operator. Q maps u to Q, a function space over D. Since 

minefw) = min min e(u), (7) 

ueu) U&fl Qu=U 

our macroscopic variational problem is given by 



where 



min E(U) 
t/efi 



E(U) = mm e(u) (9) 

Qu=U 



For dynamic problems, let us denote by s(t), the evolution operator for the microscopic 
process. In general {s(t),t > 0} forms a semi-group of operators. This semi-group may 
be generated by a set of differential equations, a Markov process, or a discrete dynamical 
system. There are at two important time scales in our problem: tn, the relaxation time 
scale of the microscopic process, and £m, the macroscopic time scale of interest. Our basic 
assumption that there exists a well-defined macroscopic model over the macroscopic time 
scale of interest implies that tn <C tjvf. 

Denote by R an appropriately chosen reconstruction operator. Given U & Q, let 

S(t)U = Qs(t)RU (10) 
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Obviously S(t)U depends on R as it is denned. However, since £r <C £m, we expect that 
the dependence on R diminishes for t ^> Therefore we can define the macroscopic model 
approximately as 

U t = F(U) (11) 

where 

v ' At 
with appropriately chosen At, t^ <^ At <^ t M . 



4 The Structure of HMM 

Our basic numerical strategy is now as follows. We will work with a macroscopic grid that 
adequately resolve the macroscopic problem, but do not necessarily resolve the microscopic 
problem, and we will attempt to solve directly the macroscopic model (|ll]) and (|8]). 

There are two main components in the heterogeneous multiscale method: An overall 
macroscopic scheme for U and estimating the missing macroscopic data from the microscopic 
model. 



4.1 The Overall Macroscopic Scheme 

The right overall macroscopic scheme depends on the nature of the problem and typically 
there are more than one choice. For variational problems, we can use the standard finite 
element method. In fact our examples in the next section use the standard piecewise linear 
finite element method. For dynamic problems that are conservative, we may use the methods 



developed for nonlinear conservation laws (see, e.g. [|44|). Examples include the Godunov 
scheme, Lax-Friedrichs scheme, and the discontinuous Galerkin method. For dynamic prob- 
lems that are non-conservative, one could simply use a standard ODE solver, such as the 
forward Euler or the Runge-Kutta method, coupled with the force estimator that we discuss 
below. 



4.2 Estimation of the Macroscopic Data 

After selecting the overall macroscopic scheme, we face the difficulty that not all data for the 
macro scheme are available since the underlying macro model is not explicitly known. The 
next component of HMM is to estimate such missing data from the microscopic model. This 
is done by solving the micro model locally subject to the constraint that Qu = U where Q is 
the approximation to Q and U is the current macro state. For example, for the variational 
homogenization problem, the missing data is the stiffness matrix for the macro model. As we 
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explain in the next section, this data can be estimated by solving the original microscopic 
variational homogenization problem on a unit cell in each element of the triangulation, 
subject to the constraint that Qu = U . For dynamic problems, such data can be estimated 
from a Godunov procedure, namely, that we first reconstruct the micro state from U, and 
evolve the micro state subject to the constraint that Qu = U, and then estimate the missing 
data from u. The missing data can be either the forces or fluxes, or a part of the forces or 
fluxes such as the eddy viscosity term in a turbulence model. We also have the option of 
carrying out a number of such microscopic calculations (e.g. with different reconstruction 
or different realization of the randomness) and extract a more accurate estimate from the 
collection of microscopic calculations. 



4.3 Examples 

To illustrate the selection of the macroscale scheme and the estimation of missing macroscale 
data from microscale models, we will discuss some examples in more detail. 

1. Variational Problems 

Examples include 



1. 



where the multiscale nature of the problem is contained in the tensor a £ (x, u) = 
(afj(x, u)) which can be of the form 

(a) a £ (x, u) = a (x, |), where a(x, y) is smooth and periodic in y with period [0, l] d . 
This is the classical homogenization problem we discussed earlier. 

(b) a £ (x, u) = a (x, |), where a(x, y) is random and stationary in y. This can be used 
to model random medium. 

(c) a £ (x,u) = a(x,u, |), where a(x,u,y) is smooth. The dependence on u makes 
this problem nonlinear. The dependence on y can be either periodic or random 
stationary. 

The macroscale problem is of the type 



mm 

ugh. 



in f \-A(x,U,VU) - f(x)U(x)\ dx 
1{d)Jd 12 J 



2. Atomistic models of crystalline solids: 

min > V(xi — Xj) 
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subject to loading conditions, where V is a pairwise atomistic potential, Xi = yi + Ui,yi 
is the position of the i-th atom before deformation, Ui is the displacement of the z-th 
atom. The macroscale problem is of the type considered in nonlinear elasticity 

min / /(W) 

u J D 

where U is the macroscale displacement field. 

For these problems, we can choose the macroscale scheme to be the standard finite element 
method over a macroscale triangulation. The macroscale data that we need to estimate is 
either J D A(x, U, VU)dx or J D f(VU)dx for U G Vh, the finite element space. These can be 
approximated via the following steps. 

1. For each element A, approximate f.A(x,U,'VU)dx or f.(VU)dx by a quadrature 
formula. 

2. For each quadrature nodes X{ 6 A, approximate A(x,U,VU)(xi) or f(VU)(xi) by 
minimizing the original microscale problem over a micro-cell A^, subject to the con- 
straint that J A u(x)dx — J A U(x)dx, f A Vu(x)dx = J A VU(x)dx, with appropri- 
ate changes for the atomistic problem. For the periodic homogenization and crystalline 
solids problems, A^ can be chosen to be a unit cell around Xi, if we replace the con- 
straint by a periodic boundary condition or the Cauchy-Born rule, as we explain in 
the next section. For the stochastic homogenization problem, A Xi should be larger 
than the correlation length. In this case, it may also be of advantageous to perform 
ensemble averages over several realizations of a [x, |V 

2. Dynamic Problems of Conservative Type 

Examples include 

1. 

d t u e = V • (a £ (x,u)Vu £ ) 
where {a e (x,u)} is as discussed above. 

2. 

d t u £ + V ■ {a £ (x)u) = 
where a e (x) = a (x, |) , a(x, y) can either be periodic or stochastic stationary in y. 

3. Kinetic models such as the Boltzmann or BGK equations. 

4. Molecular dynamics of the type discussed in Section 2. 

5. Spin-exchange models via Kawasaki dynamics |64 . 
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Other examples may include models of phase segregation, mixtures of binary fluids, elastic 
effects, etc. The macroscale models are of the type 



U t + V • J = (12) 

where U is in general a vectorial macroscale variable, J may depend on x, U, VU, etc. 

The macroscale scheme can be either a finite volume method, such as the Godunov 
scheme, or a finite element method, such as the discontinuous Galerkin method. We will 
discuss here the finite volume method. HMM based on the discontinuous Galerkin method 



is considered in 23 



The missing macroscale data for a finite volume method for (12) is the macroscale flux J 
at the cell boundaries denoted by { J. + i }. They can be estimated by the following "Godunov- 
like" procedure 

1. Select a microcell A^ + i around the cell boundary at Xj + i- 

2. From {U™}, reconstruct the microstates {u} on {A J+ i}. u should be consistent with 
{[/"} in the sense that Qu = U n , where Q is the approximation of Q restricted to 

3. Evolve the microstate u(t) using the microscale model inside {A J+ i }, with initial state 
{u}, and subject to the constraint that 

Quit) = U 



4. Evaluate the macroscale flux {J J+ i} using {u(t)}. 

The constraint Qu = U requires some additional comment. Take the example of molecular 
dynamics. If we would like to capture the macroscale behavior at the level of Euler's equa- 
tions, the constraint is simply that the average mass, momentum and energy should be given 
by the prescribed macroscale values given by {U n }. If we would like to capture the viscous 
or higher order effects, we also need to constrain the system such that the average density, 
momentum and energy gradients be given by the macroscale values. This is less convenient 
to implement, particularly if higher order gradients are required. The discontinuous Galerkin 
formulation proposed in [p3j avoids this difficulty 

The rules for selecting {A- + i} is the same as for the variational problems. As usual for 
periodic homogenization and crystalline solids problems, A - + i can be chosen to be the unit 
cell. 

3. Dynamic Problems of Nonconservative Type 

Examples include 
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d t u £ = J2 a lj( x > u ) ° " 



dxidxj 



where {a e (x,u)} is as discussed before. 



2. Spin flip models |£| that leads to Ginzburg-Landau type of equations. 
In this case, we write the macroscale model as 

U t = F{U) 

where F(U) can be a nonlinear operator acting on U. For the macroscale scheme, we choose 
an ODE solver on a grid, such as forward Euler or Runge-Kutta, and we need to estimate 
F{U) on the macro grid. 

For each macro grid point Xj, we again select a microcell Aj around Xj. The principle for 
selecting Aj is the same as before. From {Uj 1 }, we construct a piecewise polynomial of k-th 
order in Aj denoted by U™(x). The rest of the steps are the same as that for the conservative 
systems. We note that the constraint Qu = U can be interpreted as 

{u(x) - m(x))x m dx = 

for < m < k. 

4. Macroscale Markov Chains 

When the macroscale process is a Markov chain, it is natural to use a kinetic Monte Carlo 
method as the macroscale scheme. The missing data might be the transition rates between 
macro states. Estimating such data is a rather non-trivial task. It is discussed in |26| . 



In the following, for dynamic problems we will concentrate on the simplest case when 
the macroscopic scheme itself is a Godunov scheme and the missing data is the macroscopic 



forces. Extension to more general situations will be studied in [22 



5 Compression Techniques 

The key numerical problem now is how to construct approximations to E(U) and F(U). In 
this section, we review numerical techniques for efficiently approximating E(U) and F{U) 
by exploiting the separation of spatial/temporal scales. 

5.1 Compression in the Spatial Domain 

If the macroscopic and microscopic spatial scales are separated, we can effectively reduce 
the approximation of E(U) and F(U) to a unit of microscopic size on each macroscopic cell. 
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Such an idea is embodied in e.g. the quasi-continuum method. We will illustrate this with 
some examples. 

Example 1. The Variational Homogenization Problem 

Consider the variational problem 

2 E J D a - S) dx-^ - J D (13) 



where as usual we assume that a(x, y) is smooth and periodic in y with period I = [0, l] d . 
Let Th be a macroscopic triangulation of D and Vh G Hq(D) be the standard piecewise 
linear finite element space over Th- For u G Hq(D) = lu, define Qu — U G V# = f2, if 



Vitcte = / Web (14) 

K JK 

for all triangles G T#. -E^t/) as defined in (§) involves nonlocal coupling of all the 
triangles. However, we can approximate E(U) efficiently if e <C 1. This is done as follows. 
Given U £ Vh and K £ T H , denote by xk the center of mass of K, and w# the solution of 
the problem 



^■'x K +eI ,• „• ^ £ ' OXi OXj 



subject to the condition 

u(x) — U(x) is periodic with period el. (16) 

Let 

^-^jL-S^SS* (17) 

where \K\ is the volume of K, we then approximate S(C7) by 

EM = IY1 Ak ^ U ) - I U(x)f(x)dx (18) 

K ^ D 

In this example, the computation on the microscale model is reduced to a microscopic unit- 
cell problem on each macroscopic element. 

The complexity of this method is comparable to solving directly the homogenized equa- 
tion by evaluating the coefficients of the homogenized equations on each element. This is 
the minimal one can hope for. However, our method differs from solving the homogenized 
equation in one essential aspect: Our method is based on the finite e-microscale model, not 
the homogenized equation which represents the e —>■ limit. Consequently our method can 
be readily extended to more complex problem such as the nonlinear homogenization problem 
at essentially the same cost. 
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Example 2. Quasi- Continuum Method |6z 

This is a way of doing nonlinear elasticity calculations using only atomic potentials. 
Denote by xi, x 2 , . . . xjy the positions of all the individual atoms in a crystal. At zero tem- 
puerature, the position of the atoms are determined by 



N 



mm 



V(x 1 ,x 2 , • • .x N ) - ^2f(xi)ui > (19) 



i=l 



subject to appropriate boundary conditions. Here / is the external force, Ui is the displace- 
ment of the z-th atom, V is the interaction potential between the atoms. 

Quasi-continuum method works with a macroscopic triangulation of the crystal and a 
standard piecewise linear finite element space for the displacement field. The compression 
opeartor is defined in a similar way as in ([14]): Q u = U G Vjj if 



< sju > K = — I VUdx (20) 
\K\ J K 



for all K G Th, where < Vw >k denotes the average strain of the atoms on the element K. 
Having defined Q, E(U) is defined as in Section 3. 

To approximate E(U), Tadmor et.al. uses the Cauchy-Born rule. Given U G Vh, let 
ex{U) be the potential energy of a unit cell of the crystal subject to the uniform strain VZ7 
on K. Let 

E(U)= J2n K e K (U)- [ f(x)U(x)dx (21) 
k Jd 

where % is the number of unit cells on K. E(U) is the approximation of E(U). 

Quasi-continuum method contains an additional element for dealing with defects and 
interfaces in crystals, i.e. type A problems, by replacing the Cauchy-Born rule with a full- 
atom calculation on elements near the defects and interfaces. We will return to this later. 



5.2 Compression in the temporal domain 

By resorting to the microscopic model in order to simulate the macroscopic dynamics, we are 
forced to resolve the microscopic times scales which are not of interest. This is particularly 
expensive if £r tu- However in this case we can explore this time scale separation to 
reduce the computational cost in the temporal domain. 

It is helpful to distinguish two different scenarios by which relaxation to local equilibrium 
takes place. For some problems, such as the parabolic homogenization problem (^3) and the 
Boltzmann equation, we have strong convergence to equilibrium. No temporal or ensemble 
averaging is necessary for the convergence of the physical observables. For other problems, 
such as the advection homogenization problem and molecular dynamics, convergence to 



13 



equilibrium is in the sense of distributions, i.e. physical observables converge to their local 
equilibrium values after time or ensemble averaging. Let us express the approximate F(U), 
called a F-estimator, in the form 



F(U)=Qj2^f« 



where the weights {4>j} should satisfy 



Eft 



Uj is the computed microscopic state at microscopic time step j, and u = RU where R 
is some reconstruction operator. The selection of the weights in the F-estimator crucially 
depends on the nature of this convergence. In particular we note two special choices. The 
first is: ip^ = 1 and ipj = for j < k. This is suitable when we have strong convergence to 
equilibrium. The second choice is: ifij = t, for 1 < j < k. This is more suited for the case 
when we have weak convergence to local equilibrium. More accurate choices of the weights 



are discussed in |[22|| . The time interval on which the microscopic model has to be solved 
depends on how fast the transient introduced by the reconstruction step dies out. 
Consider the parabolic homogenization problem 

u \ = V • (a (x, ~) Vu £ ) (22) 

on D, with Dirichlet boundary condition u £ \qd = 0. To approximate the macroscopic behav- 
ior of u £ , we will work with a macroscopic grid of size (Ax, At). Let U = Qu e be the moving 
cell averages of u e over a cell of size Ax. Let R be the piecewise linear reconstruction. In 
one-dimension, this is RU(x) = Uj + ^^"^ ( x ~ x j), f° r x £ [jAx, (j + I) Ax). With this 
reconstruction, we proceed with the microscopic solver. Asymptotic analysis suggests that 
the relaxation time for this problem is 0(e 2 ) JJ. We plot in Figure 1 a typical behavior of 
the microscopic flux j £ (x,t) = a (x, |) \7u £ (x,t) at a cell boundary over the time interval 
[t n , t n + At] as a function of the micro time steps. It is quite clear that j £ (x, t) quickly settles 
down (after about 35 micro time steps) to a quasi-stationary value after a rapid transient. 
We obtain an efficient numerical scheme if we select this value as the macroscopic flux and 
use that to evolve U over a much larger time step At. 
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Figure 1. Computed flux r e (x,t) = a (x, I) Vu £ (x,t) clS Si function of the micro time step over 
one typical macro time step, for the parabolic homogenization a (x, |) = 2 + sin27r|. The bottom 
figure is a detailed view of the top figure for small time steps. Notice that j £ (x,t) quickly settles 
down (after about 35 micro time steps) to a quasi-stationary value after a rapid transient. 
Our next example is the advection homogenization problem 

u\ + V • (a (x, ^ u e ) = (23) 

in one-dimension. We assume a(x, y) > a > 0. We proceed as before, except that we 
take a piecewise constant reconstruction. In contrast to the previous example, the temporal 
oscillations in the solutions of (^) do not die out. This is reflected in Figure 2 where we 
plot the microscopic flux j £ (x,t) = a (x, -j u £ (x,t) over the time interval [t n ,t n + At] as a 
function of the microscale time steps. j £ remains oscillatory throughout the time interval. 
Nevertheless, if we plot the time average 

i r tn+t ( t\ 

j(x, t) = - J K[l--) f{x, r)dr, K(t) = 1 - cos 2nr (24) 

as shown in Figure 3, we see that it settles down to a quasi-stationary value on a time scale 
of O(e). 
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Figure 2. Top figure: Computed flux j £ (x,t) = a (x, |) u e as a fraction of the micro time step 
over one macro time step for the convection homogenization problem (p3[). Bottom figure: Time 
averaged flux j(x,t) as a function of the micro time step. 

The fact that the microscopic process only has to be evolved on time scales comparable 
to tn leads to other possibilities of state space compression by neglecting the part of the 
state space which does not contribute significantly to the F-estimator. 

In summary, we can express the F-estimator at time t as 

F e {U, t) = Q A t{f(u(r)),t <r<t + At, u{t) = RU} (25) 

where u(t) is the solution of the compressed microscopic model (possibly over a truncated 
computational domain) with initial data u(t) = RU, Q A t is the numerical approximation of 
the compression operator. Typically Q A t has the form 

Q At = QeQxQt (26) 

where Q e , Q x , Qt denote the compression operators over the probability, spatial and temporal 
spaces respectively. Having F £ (U,t), the macroscopic state variables can be updated using 
standard ODE solvers. The simplest example of forward Euler scheme gives 

jjn+l = Jjn + AtF £ (t/ n , t n ). (27) 

6 Stability and Accuracy of HMM 
6.1 Variational Problems 

The analysis of HMM proceeds in the same way as the analysis of traditional numerical 
methods, except we have to deal with in addition the effect of compression. For variational 
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problems, compression gives rise to additional error in the evaluation of the macroscale 
energy functional, or for linear problems, the stiffness matrix. 

Take the example of the variational homogenization problem. The main error in the 
evaluation of the stiffness matrix comes from the inconsistency at the surface of the element 
where it meets another element. This error is of the order ^-||VC/||| 2 - Consequently, one 
has 

Theorem 1. Assume that the finite element triangulation is quasi- regular, then 

\\U-Qu e \\ HHD) <c(H+^ 
where U is the numerical solution of HMM, H is the size of the macroscale element. 



6.2 Dynamical Problems 

For dynamic problems, it is helpful to define an auxiliary macroscale scheme, called the 
Generalized Godunov Scheme (GGS) in [21 . Roughly speaking, GGS is obtained if in HMM 
we replace the microscale solver in the data estimation step by the macroscale solver. Since 
the macroscale model is not explicitly known, the GGS is not a practical tool but only an 
analytical tool that is helpful for analyzing HMM. 

For example, for the parabolic and advection homogenization problems discussed earlier, 
GGS is simply the Godunov scheme on the macroscale problem with appropriate reconstruc- 
tion and approximate Riemann solvers. 

Assuming that the macroscale model is in the form of a differential equation \J t = F(U), 
we can write a one-step HMM in the form 

= U? + AtFj{U n ) 



and GGS in the form 



m +1 = m + /\tF 3 {u r - 



The basic stability result proved in is that if GGS is stable, then the HMM is stable and 
\\U n -U n \\ < C(\\U°-U-°\\ + max \\F(U k ) - F(U k ) ||) 



o< fc <si 



for nAt < T. 

The notion of stability for the GGS has to be quantified appropriately for nonlinear 
problems. See [^] for details. 
Noting that 

\\U n - Qu £ \\ < \\U n - U n \\ + \\U n - Qu £ \\ 
we now conclude that the stability and accuracy of HMM depends on 
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1. Consistency of GGS with the macroscale model. 



2. Stability of GGS. 

3. The compression error \\F(U) — F(U)\\ 
We discuss each of these in some more detail. 

1. Consistency of GGS and the macroscale model might be lost if the overall macroscale 
scheme does not probe the macroscale properties to the right level of accuracy. For 
example, if the macroscale model is hydrodynamics including viscous effect, and in 
HMM we have only probed the flux in the convective terms by using a piecewise 
constant reconstruction near the cell boundaries, neglecting the dissipative terms. This 
results in inconsistency with the macroscale model. Other such examples are discussed 
in 0. 

2. Stability of GGS usually results in the standard constraint on macro time step size. It 
may also impose constraints on the reconstruction operator. 

3. The compression error has also several sources, e.g. compressions in the temporal 
or spatial domains. The nature of the temporal compression error depends on the 
nature of the relaxation to local equilibrium of the microscale process. In the case of 
strong relaxation, no temporal averaging is necessary for macroscale data estimation. 
In the case of weak relaxation, the temporal compression error depends strongly on 
the temporal and/or ensemble averaging operator used. 

[ PT| also pointed out the importance of averaging out spatial small scales for HMM based 
on the flux-formulation. 

7 Conclusion 

There are two important questions that have to addressed in order to design efficient numer- 
ical methods that couple the macro and microscale models: 

1. What is the best way to set up the individual microscale problems? 

2. How do we couple the microscale problems together in order to simulate the macroscale 
behavior? 
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The second question is now fairly adequately addressed by HMM. The first question is 
tied more with specific applications. We have discussed a few examples. But much more 
work needs to be done in order to understand the issue in the general case. 
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